Copyright (c) 2018 Geosoft Inc.
https://github.com/GeosoftInc/gxpy
A grid is a form of spatial data that represents information (such as a gravity intensity, a magnetic reading, or a colour) at points organized as a 2-dimensional array on a right-handed cartesian plane:
The plane on which the grid is located can be oriented in three dimensions relative to a Coordinate System on the Earth. The most common grids are located on a horizontal surface relative to the Coordinate System. For example, common surfaces might be sea-level, or the ground surface, or a constant elevation above or below the ground surface, or a constant elevation. Vertical cross-sections through the Earth are oriented to be orthogonal to the surface of the Earth.
A common term used with grids is the concept of a 'grid cell'. In Geosoft's usage, grids are an array of points at a point location, and a 'grid cell' is the rectangular area that extends half-way to the neighboring grid points.
Spatial reference angles throughout the geosoft.gxpy module will consistently use an angle in degrees azimuth, which is a clockwise-positive angle relative to a coordinate system frame (North, or positive Y). The geosoft.gxapi.GXIMG class specifies the rotation angle in degrees counterclockwise-positive, and other references within the geosoft.gxapi may differ.
Grids are stored as files on the file system, and there are many common grid file formats in existence. When working with grid files in GX Developer you define the grid file format with the use of a decorator string appended to the grid file name. Geosoft supports the 11 formats (and their many derivatives) described in the Grid File Name Decorations section of the GX Developer documentation. For example:
| grid file string | grid file type |
|---|---|
| 'c:/project/mag.grd(GRD)' | Geosoft format grid. |
| 'c:/project/mag.tif(TIF)' | GeoTIF |
| 'c:/project/image.jpg(IMG;T=5)' | jpeg image file |
| 'c:/project/mag.grd(GRD;TYPE=COLOR)' | Geosoft colour grid |
See also: Tutorial Page
In [1]:
import geosoft.gxpy.gx as gx
import geosoft.gxpy.grid as gxgrid
import geosoft.gxpy.utility as gxu
from IPython.display import Image
gxc = gx.GXpy()
url = 'https://github.com/GeosoftInc/gxpy/raw/9.3/examples/tutorial/Grids%20and%20Images/'
gxu.url_retrieve(url + 'elevation_surfer.GRD')
Out[1]:
We will start with a common simple task, converting a grid from one format to another. Geosoft supports many common geospatial grid formats which can all be openned as a geosoft.gxpy.grid.Grid instance. Different formats and characteristics are specified using grid decorations, which are appended to the grid file name. See Grid File Name Decorations for all supported grid and image types and how to decorate the grid file name.
Problem: You have a grid in a Geosoft-supported format, and you need the grid in some other format to use in a different application.
Grid: elevation_surfer.grd, which is a Surfer v7 format grid file.
Approach:
(SRF;VER=V7).gxgrid.Grid.copy class method to create an ER Mapper grid, which will have decoration (ERM).
In [2]:
# open surfer grid
with gxgrid.Grid.open('elevation_surfer.grd(SRF;VER=V7)') as grid_surfer:
# copy the grid to an ER Mapper format grid file
with gxgrid.Grid.copy(grid_surfer, 'elevation.ers(ERM)', overwrite=True) as grid_erm:
print('file:', grid_erm.file_name,
'\ndecorated:', grid_erm.file_name_decorated)
You work with a grid using a geosoft.gxpy.grid.Grid instance, which is a spatial dataset sub-class of a geosoft.gxpy.geometry.Geometry. In Geosoft, all spatial objects are sub-classed from the Geometry class, and all Geometry instances have a coordinate system and spatial extents. Other spatial datasets include Geosoft databases (gdb files), voxels (geosoft_voxel files), surfaces (geosoft_surface files), 2d views, which are contained in Geosoft map files, and 3d views which can be contained in a Geosoft map file or a geosoft_3dv file.
Dataset instances will usually be associated with a file on your computer and, like Python files, you should open and work with datasets using the python with statement, which ensures that the instance and associated resources are freed after the with statement looses context.
For example, the following shows two identical ways work with a grid instance, though the with is prefered:
In [3]:
# open surfer grid, then set to None to free resources
grid_surfer = gxgrid.Grid.open('elevation_surfer.grd(SRF;VER=V7)')
print(grid_surfer.name)
grid_surfer = None
# open surfer grid using with
with gxgrid.Grid.open('elevation_surfer.grd(SRF;VER=V7)') as grid_surfer:
print(grid_surfer.name)
One often needs to see what a grid looks like, and this is accomplished by displaying the grid as an image in which the colours represent data ranges. A simple way to do this is to create a grid image file as a png file ising the image_file() method.
In this example we create a shaded image with default colouring, and we create a 500 pixel-wide image:
In [4]:
image_file = gxgrid.Grid.open('elevation_surfer.grd(SRF;VER=V7)').image_file(shade=True, pix_width=500)
Image(image_file)
Out[4]:
A nicer image might include a neat-line outline, colour legend, scale bar and title. The gxgrid.figure_map() function will create a figure-style map, which can be saved to an image file using the image_file() method of the map instance.
In [5]:
image_file = gxgrid.figure_map('elevation_surfer.grd(SRF;VER=V7)', title='Elevation').image_file(pix_width=800)
Image(image_file)
Out[5]:
In Geosoft all spatial data should have a defined coordinate system which allows data to be located on the Earth. This also takes advantage of Geosoft's ability to reproject data as required. However, in this example the Surfer grid does not store the coordinate system information, but we know that the grid uses projection 'UTM zone 54S' on datum 'GDA94'. Let's modify this script to set the coordinate system, which will be saved as part of the ER Mapper grid, which does have the ability to store the coordinate system description.
In Geosoft, well-known coordinate systems like this can be described using the form 'GDA94 / UTM zone 54S', which conforms to the SEG Grid Exchange Format standard for describing coordinate systems. You only need to set the coordinate_system property of the grid_surfer instance.
In [6]:
# define the coordinate system of the Surfer grid
with gxgrid.Grid.open('elevation_surfer.grd(SRF;VER=V7)') as grid_surfer:
grid_surfer.coordinate_system = 'GDA94 / UTM zone 54S'
# copy the grid to an ER Mapper format grid file and the coordinate system is transferred
with gxgrid.Grid.copy(grid_surfer, 'elevation.ers(ERM)', overwrite=True) as grid_erm:
print(str(grid_erm.coordinate_system))
Coordinate systems also contain the full coordinate system parameter information, from which you can construct coordinate systems in other applications.
In [7]:
with gxgrid.Grid.open('elevation.ers(ERM)') as grid_erm:
print('Grid Exchange Format coordinate system:\n', grid_erm.coordinate_system.gxf)
In [8]:
with gxgrid.Grid.open('elevation.ers(ERM)') as grid_erm:
print('ESRI WKT format:\n', grid_erm.coordinate_system.esri_wkt)
In [9]:
with gxgrid.Grid.open('elevation.ers(ERM)') as grid_erm:
print('JSON format:\n', grid_erm.coordinate_system.json)
In [10]:
# show the grid as an image
Image(gxgrid.figure_map('elevation.ers(ERM)', features=('NEATLINE', 'SCALE', 'LEGEND', 'ANNOT_LL')).image_file(pix_width=800))
Out[10]:
In this exercise we will work with the data stored in a grid. One common need is to determine some basic statistical information about the grid data, such as the minimum, maximum, mean and standard deviation. This exercise will work with the grid data a number of ways that demonstrate some useful patterns.
The smallest code and most efficient approach is to read the grid into a numpy array and then use the optimized numpy methods to determine statistics. This has the benefit of speed and simplicity at the expense memory, which may be a concern for very large grids, though on modern 64-bit computers with most grids this would be the approach of choice.
In [11]:
import numpy as np
# open the grid, using the with construct ensures resources are released
with gxgrid.Grid.open('elevation_surfer.grd(SRF;VER=V7)') as grid:
# get the data in a numpy array
data_values = grid.xyzv()[:, :, 3]
# print statistical properties
print('minimum: ', np.nanmin(data_values))
print('maximum: ', np.nanmax(data_values))
print('mean: ', np.nanmean(data_values))
print('standard deviation: ', np.nanstd(data_values))
Many Geosoft methods will work with a geosoft.gxpy.vv.GXvv, which wraps the geosoft.gxapi.GXVV class that deals with very long single-value vectors. The Geosoft GXVV methods works with Geosoft data types and, like numpy, is optimized to take advantage of multi-core processors to improve performance. The pattern in this exercise reads a grid one grid row at a time, returning a GXvv instance and accumulate statistics in an instance of the geosoft.gxapi.GXST class.
In [12]:
import geosoft.gxapi as gxapi
# the GXST class requires a desktop license
if gxc.entitled:
# create a gxapi.GXST instance to accumulate statistics
stats = gxapi.GXST.create()
# open the grid
with gxgrid.Grid.open('elevation_surfer.grd(SRF;VER=V7)') as grid:
# add data from each row to the stats instance
for row in range(grid.ny):
stats.data_vv(grid.read_row(row).gxvv)
# print statistical properties
print('minimum: ', stats.get_info(gxapi.ST_MIN))
print('maximum: ', stats.get_info(gxapi.ST_MAX))
print('mean: ', stats.get_info(gxapi.ST_MEAN))
print('standard deviation: ', stats.get_info(gxapi.ST_STDDEV))
A grid instance also behaves as an iterator that works through the grid data points by row, then by column, each iteration returning the (x, y, z, grid_value). In this example we will iterate through all points in the grid and accumulate the statistics a point at a time. This is the least-efficient way to work through a grid, but the pattern can be useful to deal with a very simple need. For example, any Geosoft supported grid can be easily converted to an ASCII file that has lists the (x, y, z, grid_value) for all points in a grid.
In [13]:
# the GXST class requires a desktop license
if gxc.entitled:
# create a gxapi.GXST instance to accumulate statistics
stats = gxapi.GXST.create()
# add each data to stats point-by-point (slow, better to use numpy or vector approach)
number_of_dummies = 0
with gxgrid.Grid.open('elevation_surfer.grd(SRF;VER=V7)') as grid:
for x, y, z, v in grid:
if v is None:
number_of_dummies += 1
else:
stats.data(v)
total_points = grid.nx * grid.ny
# print statistical properties
print('minimum: ', stats.get_info(gxapi.ST_MIN))
print('maximum: ', stats.get_info(gxapi.ST_MAX))
print('mean: ', stats.get_info(gxapi.ST_MEAN))
print('standard deviation: ', stats.get_info(gxapi.ST_STDDEV))
print('number of dummies: ', number_of_dummies)
print('number of valid data points: ', total_points - number_of_dummies)